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Particle stabilized emulsions are ubiquitous in the food and cosmetics industry, but our un- 
derstanding of the influence of microscopic fluid-particle and particle-particle interactions on the 
macroscopic rheology is still limited. In this paper we present a simulation algorithm based on a 
multicomponent lattice Boltzmann model to describe the solvents combined with a molecular dy- 
namics solver for the description of the solved particles. It is shown that the model allows a wide 
variation of fluid properties and arbitrary contact angles on the particle surfaces. We demonstrate 
its applicability by studying the transition from a "bicontinuous interfacially jammed emulsion gel" 
(bijel) to a "Pickering emulsion" in dependence on the contact angle, the particle concentration, 
and the ratio of the solvents. 

PACS numbers: 47.11.-j 47.55.Kf, 77.84.Nh, 



I. INTRODUCTION 

Using particles in a manner similar to surfactants in or- 
der to stabilize emulsions is very attractive in particular 
for the food-, cosmetics-, and medical industry to stabi- 
lize, e.g. barbecue sauces and sun cremes or in order to 
produce sophisticated ways to deliver drugs at the right 
position in the human body. The microscopic processes 
leading to the commercial interest can be understood by 
assuming an oil-water mixture. Without any additives, 
phase separation would take place and the oil would float 
on top of the water. Adding small particles, however, 
causes these particles to diffuse to the interface which is 
being stabilized due to a reduced surface energy. If for 
example individual droplets of one phase are covered by 
particles, such systems are also referred to as "Picker- 
ing emulsions" and are known since the beginning of the 
last century [l], Q . Particularly interesting properties of 
such emulsions are the blocking of Ostwald ripening and 
the rheological properties due to irreversible particle ad- 
sorption at interfaces as well as interface bridging due to 
particle monolayers [§-[1]. 

Recently, there has been a growing interest in particles 
suspended in multiphase or multicomponent flows 0, @- 
8], which led to the discovery of a new material type, 
the "bicontinous interfacially jammed emulsion gel" (bi- 
jel) @. The existence of the bijel was predicted in 2005 
by Stratford et al. [ljj [HI an d experimentally confirmed 
by Herzig et al. in 2007 [12]. In contrast to Picker- 
ing emulsions which consist of unconnected particle sta- 
bilised droplets distributed in a second continuous fluid 
phase, the bijel shows an interface between two continu- 
ous fluid phases which is covered by particles. 

Since the particles used for stabilization have a larger 
size than surfactant molecules and do not present any 
amphiphilic properties, concepts developed for the de- 
scription of surfactant stabilized systems are often not 
applicable. Instead, theoretical models have to be devel- 



oped and experiments have to be performed which con- 
sider the specific properties of particle-stabilized systems. 
These include the particle's contact angle, the strong 
interparticle capillary forces, or the pH value and elec- 
trolyte concentration of the solvents 03 • How- 
ever, even today our quantitative understanding of solid 
stabilized emulsions is still far from satisfactory. 

Computer simulations arc promising to understand 
the dynamic properties of particle stabilized multiphase 
flows. However, the shortcomings of traditional simula- 
tion methods quickly become obvious: a suitable simu- 
lation algorithm is not only required to deal with sim- 
ple fluid dynamics but has to be able to simulate sev- 
eral fluid species while also considering the motion of 
the particles and the fluid-particle interactions. Some 
recent approaches trying to solve these problems utilize 
the lattice Boltzmann method for the description of the 
solvents [l5[ ■ The lattice Boltzmann method can be seen 
as an alternative to conventional Navier Stokes solvers 
and is well established in the literature. It is attractive 
for the current application since a number of multiphase 
and multicomponent models exist which are comparably 
straightforward to implement [1614221 ]. In addition, the 
method has been combined with a molecular dynamics 
algorithm to simulate arbitrarily shaped particles in flow 
and is commonly used to study the behavior of particle- 
laden single phase flows |23l426| . 

A few groups combined multiphase lattice Boltzmann 
solvers with the known algorithms for suspended par- 
ticles [l(| H3|- I n this paper we follow an alternative 
approach: we present a method based on the multicom- 
ponent lattice Boltzmann model of Shan and Chen [l6| 
which allows the simulation of multiple fluid components 
with surface tension. Our model generally allows ar- 
bitrary movements and rotations of arbitrarily shaped 
hard shell particles. It does not require fluid-filled par- 
ticles and thus does not suffer from unphysical behavior 
caused by oscillations of the inner fluid |28| . Further, it 
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allows an arbitrary choice of the particle wettability - 
one of the most important parameters for the dynamics 
of multiphase suspensions [y, 0] • 

The remainder of the paper is organised as follows: 
after a description of the Shan-Chen approach for mul- 
ticomponent lattice Boltzmann simulations and an ex- 
tension of the lattice Boltzmann method to simulate sus- 
pensions, a way to combine the two methods is proposed. 
The influence of the parameters of the model on the con- 
tact angle as a measure of wettability is studied in the 
following section. Then the suitability of the new method 
is tested by performing a detailed study of the formation 
of bijcls and Pickering emulsions. 



II. THE MULTICOMPONENT LATTICE 
BOLTZMANN MODEL 

The dynamics of the fluid solvents is simulated by a 
multicomponent lattice Boltzmann model following the 
approach of Shan and Chen [l6| . Here, each component 
follows a lattice Boltzmann equation 



/F(x + c i; i+l) = /r(x,i) + 0^(x,i), 



(1) 



where (x, t) is the single particle distribution function 
for component c in the direction Ci (i = 1, . . . , N) at a 
discrete lattice position x and at timestep t. In this work 
we use exclusively the so-called D3Q19 implementation, 
where N = 19 velocities are used on a three dimensional 
lattice. For simplicity, the length of a timestep and the 
lattice constant are set to 1, i.e. all units are given in 
lattice units if not stated otherwise. is the Bhatnagar- 
Gross-Krook (BGK) collision operator [29l |. 



fi?(x,t) 



f?(x,t)-fP (p c (x,t),u c (x,i)) 



, (2) 



which is a relaxation towards the local equilibrium dis- 
tribution function 
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on a time scale given by the relaxation time r c [30l | . Here, 
p c (x, t) = po J^i ft ( x ; t) is the fluid density with ref- 
erence density po and u = u c (x, t) is the macroscopic 
bulk velocity of the fluid, given by /O c (x, i)u c (x, t) = 
(x, t)cj. d are the coefficients resulting from the 
velocity space discretization and c s = 1/ \/3 is the speed 
of sound, both of which are determined by the choice of 
the lattice. The kinematic viscosity of the fluid is given 
by ^ = c 2 (r c - 1/2). 

The interaction between fluid components c and c' 
is introduced as a self-consistently generated mean field 
force 



F c (x,i) = -vf c (x,0 



E 



* c '(x',;)(x'-x) , ( 4 ) 



where x' are the nearest neighbors and ^ c (x) is the so- 
called effective mass, which can have a general form for 
modeling various types of fluids. We choose 



* c (x, t) = po 1 - cxp - 



P c (x,t) 



Po 



(5) 



g cc i is a force coupling constant whose magnitude controls 
the strength of the interaction between components c, c! 
and is set positive to mimic repulsion. The dynamical 
effect of the force is realized in the BGK collision operator 
by adding to the velocity u in the equilibrium distribution 
the increment 



Au c (x,t) 



r c F c (x,i) 
P c (x,£) 



(6) 



The force also enters the calculation of the actual macro- 
scopic bulk velocity as [U HH 

u c (x,i) = MM^ + i r(X]i) . (7 ) 
p c (x, t) 2 

In this paper two fluids with identical properties are 
used which are called "blue" ("b") and "red" ("r"). To 
simplify statements about the fluid ratio at a certain po- 
sition an order parameter 



.(x,t) = p r (x,t)- i o b (x,t) 



(8) 



is introduced. The Shan-Chen model is a diffuse inter- 
face method, where interfaces between different fluids are 
about four lattice sites wide. For the analysis below we 
define the interface position to be located where the order 
parameter vanishes. 



III. SUSPENDED PARTICLES 

Pioneering work on the development of an extension 
to the lattice Boltzmann method to incorporate parti- 
cles has been done by Ladd et al. (23l - [25| . The method 
has been applied to suspensions of spherical and non- 
spherical particles by various authors [l(| [H, H3, 
Recently, the inclusion of Brownian motion was revis- 
ited and clarified in more detail [H, [35| . The suspended 
particles are assumed to be homogeneous spheres with 
radius r par . In our implementation of the method New- 
ton's equations for the momentum 



du r 



df 



and the angular momentum 



dt 



(9) 



(10) 



are solved with a leap frog integrator to simulate their 
behavior. Here, F is the force acting on a particle, m 
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is its mass and u par its velocity. D is the torque, J the 
moment of inertia and ui the angular velocity. 

The particles are discrctized on the lattice and interac- 
tions between the fluid and the particles are introduced 
by marking all lattice sites that are inside the particle as 
solid nodes for the fluid, at which bounce back boundary 
conditions are applied [251 ] . Bounce back boundary con- 
ditions reflect the incoming distributions at a site back 
to where they came from, so that the streaming step is 
modified to 

/f(x,t+l) = /F(x-c ij t), (11) 

for all i where x — Cj is not a solid node to 

ft(x,t+l) = ft(x,t), (12) 

for all i where x— c* is a solid node. Here, i' is defined as 
the index corresponding to c,j = — c^. This results in a 
no-slip boundary located halfway between the fluid and 
the solid node. The change of momentum of the fluid 
that is reflected at the boundary has to be compensated 
by a momentum change of the particle itself as given by 



Ap(t) = 2-p c (x,t)a 



(13) 



As we assume the length of a time step dt to be 1 this 
corresponds to a force 



and a torque 



F(t) = 2-p c (x,i)c, 



D(i) = F(t) x r(t), 



(14) 



(15) 



on the particle. Here, r(t) is the vector pointing from the 
center of the particle to the site of the reflection. Since 
the particles are not stationary but move over the lattice, 
the bounce back rule does not correctly reproduce the 
velocity of the reflected fluid. It is therefore corrected as 

ft (x, t + 1) = ft, (x, t) - \p c (x, t) u surf (x, t) ■ a> , (16) 

where u sur f (x, t) is the velocity of the particle surface on 
which the fluid is reflected. This effect also leads to a 
change in transfered momentum, so that the force acting 
on the particle is 

F(t) = \2p c (x,t) - ^p c (x,t)u smi (x,t) -c^ cv. (17) 

When the particle moves over the lattice, individual lat- 
tice sites can be either occupied by it in front of the par- 
ticle or be released at its back. In case of newly occupied 
sites, the fluid on the site is deleted and its momentum 
transfered to the particle by adding 



F(i) = -p c (x,t)u c (x,t). 



(18) 



In case of the particle vacating a lattice site new fluid is 
created with the initial fluid density and the velocity of 
the particle surface u sur f(x, t) at the corresponding site: 

ft(x, t) = tf nit • ff (u surf (x, f), p(x, t)) . (19) 



To satisfy conservation of momentum this again leads to 
a force on the particle, 



F(t) = /4i t u surf (x,i). 



(20) 



Interactions between particles can be taken into ac- 
count similar to standard molecular dynamics implemen- 
tations. In the current paper, we only consider Hertzian 
contact forces to mimic hard spheres and a lubrication 
correction to correct for the limitations of the lattice 
Boltzmann method to describe the hydrodynamics prop- 
erly on scales below the lattice resolution. When particles 
collide the resulting forces are derived from the Hertzian 
potential (3|| 



t t ,s I A^Hcrtz ■ (2r par - r) 2 for r < 2r par , . 
Wr)= else, (21) 



with Anertz being a constant [37]]. When two particles 
move towards each other the lubrication interaction be- 
tween them results in a force separating the particles. If 
there is not at least one lattice site between the particles 
to resolve the flow, this force is not properly reproduced 
by the simulation and a lubrication correction has to be 
added as given by 



'lub = 



67r^ c r^ ar r 
(2r par ) 2 M 



(ui - u 2 ) 



1 



r - 2r 



- 1 



par 



(22) 

where r is the vector connecting the particle centers and 
111,2 is their respective velocity (25l |38|. It is introduced 
at particle distances smaller than | lattice units and lim- 
ited to it's value at a distance of of a lattice unit 
to avoid numerical instabilities due to the divergence at 
|r| = 2r par . 



IV. PARTICLES IN MULTICOMPONENT 
FLUIDS 

In order to develop a simulation algorithm for par- 
ticles in multicomponent flows the previously described 
methods can be combined as described in the current sec- 
tion. First, when extending the coupling between parti- 
cles and fluid to multiple components the treatment of 
lattice sites that are uncovered by the moving particle 
has to be adapted. In the original algorithm for a single 
fluid, such sites are filled with the initial fluid density 
pf nit . However, re- initializing such sites with multiple 
fluid components can lead to artefacts since it is not a 
priori the case that the correct fluid composition should 
correspond to the initial state of the simulation. For ex- 
ample, one kind of fluid could appear in a region where 
only the other fluid is present. To prevent such artefacts 
we use the average surrounding fluid density 



p c {x,t) = 



1 



A^NP 



(23) 



4 



where inp are all indices i for which x + C; is a non- 
particle site. Nhp is the number of these sites. Similar 
to the method described in [13, , the uncovered sites 
are re-initialized with 

p£ ew (x,i)=p c (x,t) (24) 

following a similar approach as in Eqs. IT51 andl2Hl i.e. 

/P(x, t) = pLw(x, t) ■ fP (pf lcw (x, t), u surf (t)) (25) 

and 

F = ^^ew(x,t)u surf (t). (26) 

c 

The second modification of the original algorithms is 
required to correctly take into account the effect of the 
fluid-fluid interaction forces on the fluid in the direct 
vicinity of a particle. For the calculation of the forces 
between different fluid components, also the empty lat- 
tice sites inside a particle are considered if one follows 
Eq. |U Since there are no Shan-Chen forces acting in the 
direction from the particle to the fluid, the fluid forms a 
layer of increased density around the particle. To avoid 
this artefact, the outermost layer of lattice sites inside 
the particle is not kept empty, but is filled with a virtual 
fluid density which is equivalent to the average of the 
surrounding densities p: 

p c (x,t)=p c (x,i) (27) 

This virtual fluid inside the particles does not follow the 
lattice Boltzmann equation, i.e. the advection and colli- 
sion steps are not applied. 

Further, the Shan-Chen like force acting from the fluid 
surrounding a particle on the particle itself has to be 
accounted for. This is implemented by summing up all 
Shan-Chen forces 

fw = EE fc ( x .') ( 28 ) 

X c 

acting on every lattice site x inside the particle. This 
force and the corresponding torque are then added to the 
particle within the molecular dynamics algorithm. The 
forces on the fluid outside the particles are calculated as 
before with the virtual fluid being treated like a regular 
fluid in the Shan-Chen force computation. This leads to 
a balanced force on the fluid sites near the particle sur- 
face and therefore prevents the formation of a layer of 
increased density. This is demonstrated in Fig. [TJ where 
the left subfigure shows a particle with r par = 10 be- 
ing filled with virtual fluid while in the right subfigure 
the particle is not filled with a virtual fluid. The parti- 
cle is set at an interface created by two lamellae of red 
and blue fluids at the center of the shown area. Peri- 
odic boundary conditions are applied causing a second 
interface to appear at the left and right borders of the 
sketches. As we use a diffuse interface method for the flu- 
ids, the interfaces cover about four lattice sites depicted 



with correction (a) without correction (b) 




FIG. 1. ^2 c p c after 2000 timesteps in the presence of a par- 
ticle with r par = 10 at an interface at the center of the shown 
area. The cut on the left (a) shows a particle being filled with 
virtual fluid, while on the right (b) the particle is empty as 
in the original algorithm. Without the virtual fluid, the halo 
of increased density can clearly be seen, while adding the vir- 
tual fluid successfully allows to correct for this inconsistency. 
The parameters of the simulation were po = 1 and r = 1 for 
both species, the system of size 48 3 lattice sites was initially 
divided into two lamellae of width 24 lattice sites with density 
p T = p° = 0.7, respectively. All units are given in lattice units 
throughout this paper. 

by the varying grey scale. Without the virtual fluid, the 
halo of increased density can clearly be seen, while adding 
the virtual fluid successfully allows to correct for this in- 
consistency. 

The advantage of a virtual fluid inside the particles is 
that it can be utilized to modify the wettability of the 
particles. Here, we follow an approach which has been 
introduced to model hydrophobic fluid-surface interac- 
tions for studying flow in hydrophobic microchannels or 
droplets on surfaces with arbitrary contact angles [40l — 
HU . Our approach is equivalent to the method presented 
in [28[ . The Shan-Chen interaction between the parti- 
cles and the fluids can be modified by tuning the density 
of the local virtual fluids. Increasing one of them by an 
amount |Ap| causes the particle surface to "prefer" this 
fluid with respect to the other one, i.e. the repulsion be- 
tween the increased component and the unmodified one 
increases. Ap is called "particle color" and a positive 
particle color is defined as an addition of "red" fluid, i.e. 

/new = P 1 + |Ap| , (29) 

whereas a negative color corresponds to "blue" fluid be- 
ing added, i.e. 

p h aew = P h + \Ap\ . (30) 

In the next section we demonstrate that the particle color 
can be used to tune the contact angle of the particle sur- 
face at an interface in order to resemble specific fluids and 
solid materials. As an alternative to the virtual fluid, a 
modified version of Eq. 0] that takes solid lattice sites 
and fluid-surface interactions into account could be de- 
veloped, but the approach presented here is simpler to 
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implement and does not have any relevant impact on the 
performance of the code. 

The changing discretization of the particle together 
with the fluid-surface interaction force leads to slight 
mass errors during the step in which vacated lattice nodes 
are refilled with fluid. This effect is especially strong 
for small particles and high forces (large values for g cc >). 
Typical test cases have shown that the total mass after 
very long simulation times (10 7 timesteps) increases by 
about one percent. This can be explained by the simple 
interpolation for the amount of newly created fluid which 
is necessary since no analytical solution for the multicom- 
ponent Shan-Chen model is known that describes the 
density profile at interfaces. Even though the effect is 
very small, it can be suppressed if newly created fluid 
densities p° ew are scaled with a correction factor which 
depends on the total mass error Ap c up to the current 
timcstep and the total number of lattice sites in the sys- 
tem TV. This leads to a modification of Eq. [MJ 

^ew = P C (l-^0^^^). (31) 

The rate of the adaptive correction can be tuned with the 
parameter Co- Due to the very small mass error, the cor- 
rection can act very slowly, but should not be chosen too 
fast in order to avoid hysteresis effects. Further, Eq. [31] 
reduces unphysical density gradients at particle surfaces 
and thus contributes to the stability of the algorithm. 
Repeating the same test case as above with a correction 
factor of Co = 10 results in a deviation of the mass of 
0.03 percent after 10 7 timesteps. 

V. CONTACT ANGLE MEASUREMENTS 

The contact angle 9 is a common measure for the wet- 
tability of the particle by the two fluid components. The 
influence of various simulation parameters on the contact 
angle is investigated in the current section. In order to 




FIG. 2. Definition of the contact angle 9 for a particle of 
radius r par at the interface between the two fluid components. 
Ah is the distance from the particle center to the interface. 

measure the contact angle the following setup is used in 
the simulations: in z direction, one half of the lattice 



is filled with one fluid component, the other half with 
the second one. The particle is placed at the interface 
at t = and the simulation is started. The interface 
between the two components is tracked via the (linearly 
interpolated) position at which the order parameter is 
zero. Then, the contact angle 9 can be calculated by 

cos(#) = — , (32) 

where Ah is the difference between the interface position 
and the particle center in the direction perpendicular to 
the interface (cf. figure [5]). 
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FIG. 3. Contact angle versus time for a 48x48x256 system 
with particle color 0.02, r par = 10, (?br = 0.08, and fluid 
density pi n i t = 0.6. After a rapid change from about 93° to 
about 107° the contact angle stays at a constant value and 
only shows small oscillations with an amplitude of about 1°. 

To study the time-evolution of the contact angle a sys- 
tem of size 48x48x256 lattice sites is used. The particle 
has a radius of 10 lattice sites and a color of 0.02. The 
initial fluid density is 0.6 and <?br = 0.08. Figure [3] shows 
the time dependence of the contact angle. After a short, 
fast movement at the beginning of the simulation the con- 
tact angle oscillates slightly around a fixed value. Here 
and for all further graphs in this section, we average the 
contact angle over the timesteps from 6T0 5 to 9T0 5 lead- 
ing to a final value of 9 = (107.03 ± 0.26)°. The error is 
given by the standard deviation of the data. Relating the 
variation of the angle to a change of the position of the 
interface on the lattice with regards to the particle center 
results in Ah = (—2.93 ±0.04) lattice units. The error 
in the position measurement is very small with respect 
to the lattice resolution. 

Figure [4ji shows the resulting contact angle for differ- 
ent particle sizes between r par = 2 and r par = 16. For 
small particle radii the error increases substantially and 
the measured angle is not equivalent to the one measured 
for larger particles, but is up to 15° larger. For exam- 
ple, for r par = 2 the contact angle is (120.9 ± 6.0)°. For 
particles larger than r par = 5 the error stays below 1.5° 
and the measured angles arc in the range of 106 to 110°. 
Smaller particles are more susceptible to small forces be- 
cause of their smaller mass, also one has to keep in mind 
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that our lattice Boltzmann multicomponcnt model is a 
diffuse interface method. Since the interface is about 
four lattice sites wide, small particles are completely in- 
side the interface region. For particles with a non-integer 
radius the error and the angle are larger than for parti- 
cles with integer radii. This can be adhered to being a 
discretisation effect as the particles move on the lattice. 
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FIG. 4. a) Contact angle versus particle size. The system 
size is 48x48x256, the particle color 0.02, gtr = 0.08, and 
Pinit = 0.6. The average measured angle and errors given by 
the standard deviation over timesteps from 6 • 10 5 to 9 • 10 5 
are shown. The contact angle for r par = 2 is (120.9 ± 6.0)°, 
for particles larger than r par = 5 the error stays below 1.5° 
and the angles are in the range of 106° to 110° . Particles with 
a non-integer radius show larger contact angles, which can be 
adhered to a discretization effect. 

b) Contact angle versus particle color for pi n it = 0.7. The data 
can be fitted with the equation 8 — 442- Ap+90 (dashed line). 

The dependency of the contact angle 8 on the particle 
color Ap is shown in Fig. [4]b. One can see an almost linear 
relation between the contact angle and the particle color 
in the range from a color of Ap = —0.125 (contact angle 
33.2°) to Ap = 0.125 (contact angle 147.6°). Included 
in the figure is a linear fit given by 9 — 442 ■ Ap + 90 
to stress this linear behavior. The simulations with a 
particle color of Ap > 0.15 and Ap < —0.15 result in a 
detachment of the particle from the interface. Thus, it is 
possible to choose a specific particle color to obtain the 
related contact angle or to force detachment from one of 
the fluids. 

As a next step we investigate the influence of the 
strength of the fluid-fluid interaction force on the con- 
tact angle. For strong forces determined by large gbr the 
interface is well defined and the surface tension high. Low 
5br cause a low surface tension and thus a more diffuse 
interface. When the coupling constant gbr is varied, the 
contact angle 8 changes as shown in Fig. [5^,. The stronger 
the force the stronger the particle is kept at the interface. 
If the force is too weak the particle cannot be held at the 
interface anymore. For g^ r > 0.1 almost no change to the 
contact angle can be observed and 8 converges to 93.0° 
with an error smaller than 0.1°. For < 0.08 the con- 
tact angle increases dramatically until the particle does 
not stay attached to the interface at g^ Y < 0.07. On the 



one hand, a well defined contact angle and well defined 
interfaces are preferrable requiring large values of gbr- On 
the other hand, too large values can cause very high lo- 
cal flow velocities and the lattice Boltzmann method can 
become unstable. Thus, gbr should be chosen as small as 
possible. 
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FIG. 5. a) Contact angle versus g^r- For gbr > 0.1 a 
contact angle of 93.0° with an error smaller than 0.1° is 
measured. The contact angle increases from this value at 
gbr = 0.1 to 98.0° at gbr = 0.08 with an error below one de- 
gree, gbr = 0.075 results in a contact angle of (113.6 ± 0.3)° 
and for smaller gbr the particle detaches from the interface. 

b) Contact angle versus the initial fluid density pinit- 
Pinit = 0.9 results in a contact angle of 92.5°. Reducing pi n i t 
causes the contact angle to decrease to 108.4° at pi n i t = 0.55. 

As can be observed in Fig.0D, a variation of the initial 
fluid density pinit has a similar effect as a modification of 
the coupling constant on the contact angle. However, 8 
only changes from 108.5° (pinit = 0.55) to 92.5° (pinit = 
0.9), i.e. the effect is much weaker. The reason for the 
lower impact on 8 is given by our particular choice of the 
effective mass \l/ c (x) (see Eq. [5]) which causes a damping 
of the interactions for large densities. 

The knowledge of the contact angle and particle shape 
together with a measurement of the surface tension be- 
tween both fluids allows to measure the energy required 
to detach a trapped particle from an interface [f| . While 
for spherical particles it is straightforward to compute the 
detachment energy analytically, for highly anisotropic or 
complex shaped particles this is not easily possible. A 
simulation study based on the model proposed in this 
paper would allow a well founded understanding of the 
dependence of detachment energies on particle properties 
and could be compared to experimental data. 



VI. BIJEL FORMATION 

The formation of a "bijel" (bicontinuous interfacially 
jammed emulsion gel) was first predicted by Stratford et. 
al. in 2005 [Io| . As stated in the introduction, bijels can 
form when (colloidal) particles are added to a mixture of 
two immiscible fluids. During the phase separation of the 
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two fluids, the particles accumulate at the interface until 
those are fully jammed. Since the simulations performed 
by Stratford et. al. utilize a free energy based multiphase 
lattice Boltzmann model, we show in this section that 
the multicomponent model introduced in this paper is 
also able to model the formation of a bijel. We study 
the temporal development of the system and compare 
our results with the results of Stratford et al.. Further, 
we investigate the influence of the particle concentration, 
t/bn an d Pinit = YlcPinit 011 * ne bijel formation. The 
initial conditions for the simulations are as follows: an 
identical amount of the two fluid species is distributed 
randomly throughout the system. The initial positions 
of the colorless particles (0 = 90°) are also chosen at 
random. In order to keep the system size at manageable 
256 3 lattice units and to be able to simulate a significant 
number of particles, the particle radius is kept at r par = 5 
lattice units in all simulations. 

The conversion from lattice units to SI units of a sys- 
tem containing two identical fluids with the speed of 
sound (c s = 1482.35m/s) and kinematic visosity (y = 
1.004- lCT 6 m 2 /s) of water at 20° C results inatimestep of 
At = 9.14-10~ 13 s and a lattice constant of Ax = 2.35nm. 
Since we set our particle diameter to 10 lattice units, the 
physical diameter would be 23.5nm and the side length 
of a cubic simulation volume with 256 lattice units cor- 
responds to 601. 6nm. The systems presented in this sec- 
tion are simulated for 2.8 • 10 5 timesteps corresponding 
to 2.56 ■ 10~ 7 s. 
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6. Average domain size for a system with = 0.08, 
Pinit = 0.7 and a particle concentration of 20%. During the 
first 2.5 • 10 4 steps of the simulation the average domain size 
grows from below 5 to its final value of about 31 lattice units. 

As already stated correctly by Stratford et al., thermal 
fluctuations have only little effect on the phase separa- 
tion and bijel formation and can thus be ignored in the 
simulations [ToT ]. 

To analyze the development of structures in the sim- 
ulated systems, we define the averaged time dependent 
lateral domain size L(t) which consists of an average of 
its components Li along direction i = x, y, z as given by 
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is the second order moment of the three-dimensional 
structure function 
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with respect to the Cartesian component i, () denotes 
the average in Fourier space, weighted by S(k, t) and N 
is the number of nodes of the lattice, <f>' k (i) the Fourier 
transform of the fluctuations of the order parameter 
<// = (/)— ((f)) , and ki is the zth component of the wave 
vector [44|. The simulations are performed using a 256 3 
lattice, a coupling constant of gt>r = 0.08, an initial fluid 
density of pi n i t = 0.7, a particle volume ratio a of 20 
percent (about 6400 particles), a particle size of r par = 5 
lattice units, and a particle density of 1, i.e. the par- 
ticles are slightly heavier than the fluid. The time de- 
velopment of the average domain size L(t) is shown in 
Fig. El The figure clearly shows that the system comes 
to arrest after a brief period of phase separation. Dur- 
ing the first 2.5 • 10 4 timesteps the average domain size 
L(t) increases from 5 lattice units to about 31 lattice 
units and stays at that value until the end of the simula- 
tion at t = 2.8 • 10 5 timesteps. This qualitatively agrees 
to the results obtained by Stratford et al. [13]. However, 
their simulation does not converge to a fixed domain size, 
which might be caused by the thermal motion incorpo- 
rated in their model. While thermal fluctuations are not 
strong enough to detach the particles from the interface 
they might cause some local reordering of the particles 
and therefore support further domain growth. This effect 
would be favored by the small particle diameter used in 
the simulations presented in [lfj since the particle size 
is of the same order or even smaller than the interface 
thickness. 

The arrest of the phase separation process can also be 
observed by visualizing a 2D cut at z = of the order 
parameter as in Fig. [7] or a 3D visualisation of <f> as in 
Fig. [HI The differences between timesteps t — 5000 and 
t = 10000 are large while the system barely changes be- 
tween timesteps 2.5 • 10 5 and 2.8 • 10 5 . Both visualisations 
clearly demonstrate the bicontinuouity of the fluid do- 
mains. In particular Fig. [5] depicts how the particles get 
trapped at the fluid-fluid interface and cause the demix- 
ing process to stop. 

Modifying the strength of the fluid-fluid interaction 
force by varying the coupling constant gb? also influ- 
ences the resulting domain size. This is demonstrated 
in Fig. [5^, where the averaged lateral domain size L is 
shown after t = 2.8 • 10 5 timesteps and for different <7br- 
While gb r = 0.07 leads to a domain size of about 33.7 
lattice units, gbr = 0.125 results in an average size of 
28 lattice units. The differences between the spatial di- 
rections at a certain <7br are below 0.5 lattice units. A 
higher value of the coupling constant leads to stronger 
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FIG. 7. 2D cut of the order parameter <j> at z — through the 
system studied in Fig. [6] The spots where = correspond 
to the regions occupied by particles. <f> shows pronounced dif- 
ferences between t = 5000 and t = 10000, while at timesteps 
2.5 ■ 10 5 and 2.8 ■ 10 s 
served. 
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FIG. 8. (Color online) 3D visualisation of the system pre- 
sented in Figs.[6]and[7| Shown are the particles (in green/light 
gray) and the two fluids (in red/medium gray and blue/dark 



, . , , , gray, respectively). The visualizations for t = 2.5 • 10 and 

only minor rearrangements can be ob- ° •" v . ■' ' 

t = 2.8 • 10 nicely depict the bicontinuouity of the fluids and 

the attachment of the particles to the interface. 



forces attaching the particles to the interface. Therefore, 
the size of the interface increases because the particles 
cannot slightly shift away from it in order to accomo- 
date more particles on the same interfacial area. As the 
number of particles in the system is kept constant, the 
interfacial area has to increase and therefore the resulting 
domains become smaller. 

For a variation of the initial bulk density p; n it the argu- 
ments of the previous paragraph still hold. A larger value 
of pinit leads to stronger interaction forces and therefore 
to smaller structures as described above. While the cou- 
pling constant gbi directly changes the strength of the 
force pi n it affects the force only indirectly through the 
effective mass which causes the effect to be less pro- 
nounced. As depicted in Fig. [5p the average domain size 
L decreases from about 32 lattice units for pi n j t = 0.6 to 
below 30 lattice units for p lni t = 0.9. 

The connection between the area of the interface cov- 
ered by particles and the size of the resulting structures 
can be best shown by varying the particle concentration 
a. This is depicted in Fig. 1101 Increasing particle con- 
centration leads to a larger interfacial area and therefore 
to finer structures. While the average domain size of 
a system with a particle concentration of 0.15 is about 
36 lattice units this value decreases to about 22 lattice 
units for a concentration of 0.35. A too low particle con- 
centration leads to such a small stabilized surface that 
finite size effects start to appear as they arc well known 
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FIG. 9. a) Average domain size versus gbr at t = 2.8 ■ 10 5 . 
With increasing gb r L decreases from 33.7 lattice units at 
gbi = 0.07 to 28 lattice units at gtbr = 0.125. 

b) Average domain size versus pmit- pinit = 0.6 leads to a 
L=32 lattice units. A larger value of pimt = 0.9 reduces L to 
a value slightly below 30 lattice units. Error bars are given 
by the maximum deviation of L x , L y , L z from the mean. 



from lattice Boltzmann simulations of spinodal decom- 
position [44|. This can be seen for a concentration of 
5%. Here, the structure size increases drastically, also 
the average domain size is not the same for all spatial 
directions anymore and varies by about 3 lattice units. 
It is possible to fit a function of the form L = a/a + b 
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FIG. 10. Average domain size versus particle concentration 
a. Decreasing a Irom 0.35 to 0.15 leads to an increase of L 
from 21.5 lattice units to 36 lattice units. If the concentration 
is further reduced to 0.05, finite size effects start to occur. 
Also shown is L = — + 10.85, the result of a fit to the 
concentration values between 0.15 and 0.35. Error bars are 
given by the maximum deviation of L x , L y , L z from the mean. 
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t = 10000 




with a = 3.85936 and b = 10.8479 to the values where 
finite size effect do not play a role. 

VII. PICKERING EMULSIONS 
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FIG. 11. Average domain size over time for a system of size 
256 3 , <?br = 0.09, pinit = 0.66, fluid ratio 1:3, particle color 
—0.01 and particle concentration a — 0.15. After a rapid 
growth of the domain size to over 25 lattice units during the 
first 25000 timesteps the domain size increases only slowly to 
slightly over 27 lattice units at timestep 3.0 • 10 5 . 

Another well known phenomenon that can be observed 
in mixtures of immiscible fluids and particles are Pick- 
ering emulsions 0, • Here, the particles are not neces- 
sarily equally wettable by the fluids anymore. Also, the 
ratio of the amount of fluid of different species deviates 
from 1 . The result is a system where one phase is contin- 
uous while the other forms droplets which are stabilized 
by the particles. The particles prevent the droplets from 
merging when they collide and therefore stop the growth 
of the average droplet size. As before the droplet size can 



t = 250000 t = 300000 

FIG. 12. (Color online) 3D visualisation of the system de- 
scribed in Fig. 1111 The particles (green/light gray) and the 
fluids (red/medium gray and blue/dark grey) are shown. The 
particles are attached to the interface between the fluid com- 
ponents and the red fluid forms spherical droplets inside the 
continuous blue fluid. While the change from timesteps 5000 
to 10000 is significant the droplet growth is almost at rest 
between t = 2.5 • 10 5 and t = 3.0 • 10 5 . 



be measured utilizing L(t) as shown in Fig.QTJ Here, the 
lattice is 256 3 , the interaction constant <?br = 0.09, and 
the initial fluid density pi n j t = 0.66. The fluid ratio is 
1:3 and the particles with a color of A<j> = —0.01 have 
a volume concentration of 15%. As for the previously 
presented "bijels" a rapid growth of L(t) from 5 to over 
25 lattice units can be observed during the first 25000 
timesteps of the simulation. The domain growth slows 
down to a slight decelerating growth afterwards. This 
agrees qualitatively with experimental results by Arditty 
et al. Q. 

A three-dimensional visualisation of the order param- 
eter </> and the particles is shown in Fig. [T^] and accom- 
panied by a two dimensional cut of the system at z = 
in Fig. Q2] The particle covered droplets as well as the 
slowing down of droplet growth can be observed. While 
the system changes dramatically between timesteps 5000 
and 10000 almost no difference can be observed when 
comparing step 2.5 • 10 5 to step 3.0 • 10 5 . 

The influence of the interfacial area on the droplet size 
can be demonstrated by modifying the particle concen- 
tration. The resulting average domain size for different 
concentrations is shown in Fig. 1141 A higher concentra- 
tion leads to a larger stabilised interfacial area resulting 
in smaller droplets: reducing the particle concentration 
from 0.15 to 0.05 corresponds to an increase of the av- 
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FIG. 15. Average domain size after 1.5 ■ 10 5 timesteps versus 
contact angle 9. System size 256 3 , g hr = 0.08, p init = 0.7, 
fluid ratio 5:9 and particle concentration a = 0.15. Neutrally 
wetting particles lead to a small L with a large deviation 
between the spatial directions while strongly colored particles 
lead to a larger average domain size with almost no deviations. 
Error bars are given by the maximum deviation of L x , L y , L z 
from the mean. 



FIG. 13. 2D cut at z = through the system described in 
figure [TT1 at different times. Shown is the order parameter </>. 
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FIG. 14. Average domain size after 1 • 10 6 timesteps versus 
particle concentration a. System size 256 3 , ghi = 0.08, pinit = 
0.7, fluid ratio 1:3 and particle color —0.01. Increasing a leads 
to a decreasing average domain size. Also shown is a fit with 
+ 18.91. Error bars are given by the maximum 
from the mean. 



equation 
deviation of L 



L y , L z 



erage droplet size from 29 to 41 lattice units. As the 
simulated system is finite, modifying the concentration 
of particles does also change the volume of the two fluid 
components. Therefore, the inversely proportional rela- 
tion between particle concentration and droplet size as 
found by Arditty et al. 0] does not apply here, but has 
to be shifted by a constant offset. L = + 18.91 is 
found to be a good fit of the data presented in Fig. [14] 

As expected the colour and thus the contact angle of 
the particle has a drastic influence on the formed struc- 
ture. While strongly colored particles with contact an- 
gles different from 90 degrees lead to spherical droplets, 



neutrally wetting particles result in droplets that are not 
as spherical anymore. We observe structures that are 
similar to the ones found by Kim et al. [l"lT ] for their sim- 
ulation of neutrally wetting particles. These structures 
are extended in one of the directions. This results in a re- 
duction of the measured average domain size L, while the 
difference between the directions increases. This differ- 
ence is expressed through the error bar in Fig. [15] The 
values of the contact angle shown in the figure are ob- 
tained from the mapping presented in Fig. 0}d. 



VIII. TRANSITION FROM BIJEL TO 
PICKERING EMULSION 

In the current section it is demonstrated how the con- 
tact angle, the particle volume concentration and the ra- 
tio of the two fluid species determine the final state of 
the system to be a bijel or a Pickering emulsion. Phase 
diagrams depending on the various simulation parame- 
ters arc presented in Fig. [16] and [IT] In order to reduce 
the computational cost, the size of the lattice has been 
reduced to 128 3 in this section. However, by performing 
a small number of 512 3 sized simulations it has been con- 
firmed that finite size effects are still below an acceptable 
limit and do not influence the final physical state of the 
system. If the system categorizes as bijel or Pickering 
emulsion is determined visually after 3 • 10 timesteps. 
Pinit is kept fixed at 0.7 and gb r is set to 0.08 in all sim- 
ulations. 

Figure [TBI shows a phase diagram in dependence on the 
particle concentration and the contact angle (see Fig. 3Jd 
for the mapping between the particle color and the con- 
tact angle). The ratio of the two fluid species is kept 
fixed at 5:9. For contact angles larger than 90 degrees the 
system always relaxes towards a bijel, while for strongly 
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FIG. 16. System state in dependence on the contact angle 
and concentration after 3 • 10 4 timesteps. The system size 
is 128 3 , /9i n it = 0.7, and the ratio of the two fluid species is 
kept fixed at 5:9. For contact angles larger than 90 degrees 
the system always relaxes towards a bijel, while for strongly 
negative coloring a Pickering emulsion is obtained. The line 
is a guide to the eye. 

negative coloring and thus smaller contact angles a Pick- 
ering emulsion is obtained. The particle concentration, 
however, only has a minor influence on the final state. 
It can only be noticed that for small concentrations the 
formation of a Pickering emulsion is more favored for 
smaller contact angles. The line is only a guide to the 
eye since the exact position of the transition from bijel 
to Pickering emulsion would require substantially more 
data points. 

In Fig. [T7] the final system state is depicted in depen- 
dence on the contact angle and the fluid ratio. A fluid 
ratio of at least 3:4 results also for contact angles larger 
than 90 degrees in a bijel, while for a fluid ratio of 2:5 
even neutrally wetting particles arc able to stabilize a 
Pickering emulsion. As already shown in Fig. [TH] for a 
fluid ratio of 5:9 it depends on the contact angle if the 
system relaxes towards a bijel or a Pickering emulsion. 
This behavior can be explained by the interplay between 
interface curvature and interface size: it is only energeti- 
cally beneficial if the work required to maintain a curved 
interface is not larger than the cost due to the increased 
size of the interface, where the latter can be overcome by 
a higher concentration of particles at the interface. 

IX. CONCLUSION 

In this paper we proposed a new method allowing the 
simulation of particles with variable contact angle in mul- 



ticomponcnt fluid flows. We have studied the influence of 
the model parameters on the resulting fluid-particle inter- 
actions and shown that our approach is able to simulate 
the formation of "bijels" and Pickering emulsions. By 
computing phase diagrams we have demonstrated how 
the transition from bijel to Pickering emulsion is deter- 
mined by the contact angle between particle and fluids, 
1:1 3:4 5:9 2:5 
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FIG. 17. System state in dependence on the contact angle 
and the fluid ratio after 3 • 10 4 timesteps. The system size is 
128 3 , pi n it = 0.7, and the particle concentration is kept fixed 
at 0.2. For fluid ratios between 1:1 and 3:4 also for contact 
angles larger than 90 degrees a bijel is obtained. For larger 
concentration ratios it depends on the particle color or contact 
angle if a bijel or a Pickering emulsion is produced. The line 
is a guide to the eye. 



the particle concentration, and the ratio of the two fluid 
species: while the wettability of the particles and the 
fluid ratio strongly influence the transition from a bijel 
to a Pickering state, the particle volume concentration 
only has a minor impact. 
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